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We present the analytic solution of the Mode III steady-state crack in a square lattice with 
piecewise linear springs and Kelvin viscosity. We show how the results simplify in the limit of large 
width. We relate our results to a model where the continuum limit is taken only along the crack 
direction. We present results for small velocity, and for large viscosity, and discuss the structure of 
the critical bifurcation for small velocity. We compute the size of the process zone wherein standard 
continuum elasticity theory breaks down. 



I. INTRODUCTION 

The problem of the dynamics of cracks has received renewed interest recently, [1] motivated in large part by new sets 
of experiments. [2,3] These experiments have called into question some of the predictions of the traditional, continuum 
mechanics approach to fracture dynamics. The most striking experimental finding is that cracks exhibit a branching 
instability long before they reach the predicted limiting speed of advance. This instability causes increased dissipation 
and sets an effective limit on the speed of crack propagation. There are hints of such an instability in the continuum 
approach, [4] but a systematic treatment remains elusive. [5] 

One avenue of exploration that has proven fruitful is the lattice models of fracture pioneered by Slepyan [6,7] and 
further developed by Marder and collaborators. [9,10] These models, especially in the extreme brittle limit, are simple 
enough to allow comprehensive study, both analytically and by numerical simulation. The lattice models exhibit some 
novel effects, not seen in the continuum description. Foremost is the existence of arrested cracks. The lattice models 
also show instabilities at large velocities that may be relevant to the experimentally seen branching instabilities. Thus, 
it is useful to understand the lattice models in as much detail as possible. 

In a previous paper, [11] we embarked on a study of the effect of dissipation, in the form of a Kelvin viscosity, [12] on 
the behavior of steady-state cracks. We solved numerically for the dependence of velocity as a function of the driving 
displacement A. We found that dissipation acts to lower the velocity and significantly reduces the size of the lattice- 
induced small velocity unstable regime where the velocity is a decreasing function of the driving. We also showed 
that in the presence of dissipation, the stable regime is well approximated by a novel x-continuum model, wherein the 
lattice structure perpendicular to the crack is retained but along the crack is replaced by a naive continuum limit. 
We also showed that if the transverse dimension N is large, then at distances of order N the elastic fields are given 
by the results of standard continuum fracture theory. On small scales, however, there is a boundary layer where the 
discreteness of the lattice in the transverse direction is important. This boundary layer structure is all important in 
determining the velocity versus driving relation. However, as our x-continuum model demonstrated, the discreteness 
in the direction of the crack is less crucial, and primarily affects the small velocity regime. 

In this current paper, we study the large- limit of the theory. We do this first for our x-continuum model, where 
the structure of the theory is simpler. We then extend this to the full lattice model. In both cases, we present a 
^ formal Wiener-Hopf solution of the model for arbitrary N, and then take the large-A^ limit. This is in contrast to 
Q , the work of Slepyan, who, for the case of infinitesimal dissipation, solves the infinite-iV limit directly. The principal 
O ■ advantage of our method is that it allows a discussion of case of large, but finite, N. It also allows a comparison 
^ ' between the small-scale and the large-scale structure, whereas Slepyan's method only produces a solution for small 
• . to intermediate scales. Thus Slepyan must rely on an implicit matching to large scales via the stress-intensity factor, 
' as opposed to the explicit matching contained in our solution. The Slepyan method, nevertheless, by avoiding the 
neccesity of solving the finite- A'^ problem, is more easily applied to other cases, such as the mode-I problem, where 
the finite- A^ solution is not so easily obtained. 

The plan of the paper is as follows. In Section II, we describe the lattice model and the simpler x-continuum 
version. In Section III, we lay out our major results. The details of the calculation are contained in the following 
sections, first for the x-continuum problem in Section IV, and then for the lattice problem in Section V. The small 
velocity limit is studied in Section VI and the large viscosity limit in Section VII. We conclude with some comments 
in Section VIII. 
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II. DESCRIPTION OF THE MODELS 



The lattice model we study is identical to that described in our earlier work, [11]. We have a square lattice of mass 
points undergoing (scalar) displacement out of the plane. The lattice extends infinitely long in the .T-direction, with 

+ 1 rows in the y-direction. The lattice points are connected by linear "springs", with spring constant 1, to their 
nearest neighbors. The top row is displaced a fixed amount A. The bottom row is connected to a fixed line, with 
piece-wise linear springs. These springs, with spring constant fc, "crack" irreversibly if they are stretched an amoimt 
e. When k ~ 2, this model is equivalent to a system of 2N + 2 rows, loaded by ±A from top and bottom, with a 
symmetric crack running down the middle, with extension at cracking of the springs that bridge the middle being 2e. 
All the (uncracked) springs have a viscous damping r]. The equation of motion for the system is then 



for j ^ 1 with Uj,jv+i = A, and 



1 + ) ("i+ij 



(1) 



1 



V 



dt 



{ui+ij + u^-ij + Ui^2 - 3uij) - k9{e - Ui,i) ( 1 + ri— ) m^i 



(2) 



Note that in these units, the elastic wave speed is unity, so all velocities are dimensionless, expressed as fractions of 
the wave speed. 

We are interested in steady-state cracks, described by the Slepyan traveling wave ansatz. 



Uij{t) = Uj{t — i/v) 



(3) 



which implies that every mass point in a given row undergoes the same time history, translated in time. We choose 
the origin at time such that ■Ui(O) = e so that it represents the moment of cracking of the spring attached to the 
bottom row mass point. The equation of motion is best expressed in terms of the N x N coupling matrix 



MNim) = 



-(m+1) 1 



-2 1 
1 -2 



1 -2 1 
1 -2 



(4) 



The steady-state equation then reads 



d 



iijit) - e{t) ( 1 + ) A^,j'(o)%'(i) - e{-t) [i + T)-] Mj,r{k)uj,{t) 



d 



dt 



1 + ) + V^) - 2% W + Uj{t - 1/v)) = 



(5) 



We will also consider in this paper an x-continuum version of this model, where we replace the nonlocal in time 
coupling along the crack with its continuum analog 



d 



d 



Ujit) - eit) (l + V-,) Mi,j,{0)uj,{t) - 6{-t) ( 1 + ) Mj,r{k)u^.{t) 



d 



dt ' 



Uj{t) = Q 



(6) 



III. SURVEY OF RESULTS 



In this section, we survey the major results derived in the bulk of the paper. As the derivations arc exceedingly 
technical, it is useful to present the results first by themselves so that they may be appreciated without getting lost 
in a welter of technical complications. 
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Wc begin by completing the Wiener-Hopf (WH) solution of the continuous x, discrete y model, as the results are 
simpler and are a useful basis for assimilating the more complicated results of the full lattice model. The key aspect 
of the solution is the calculation of A as a function of the crack velocity v (in units where the wave speed is unity). 
We find 

which expresses A (normalized to the Griffith value 

Ag = eV2iV + 1 (8) 

at which the uncracked state becomes metastable) in terms of the wave vectors corresponding to the various normal 
models of the problem. If we label the normal mode eigenvalues of the y-coupling matrix on the uncracked side A4{k) 
by A,„, then Qi^m is the unique positive root of the dispersion relation 

TjvQ"^ + {l-v^)Q^ + {l + 7jvQ)A^^0 . (9) 

Similarly, qi^m is the unique positive root of the dispersion relation using the normal mode eigenvalue Xm of the 
cracked side A^(0). 

This formula is fairly complicated, but simplifies tremendously for the case of symmetric cracks {k = 2) in the 
macroscopic limit A'^ >> 1. Then, the product above can be performed analytically, with the simple result 



where Qi(l) is the mode associated with the highest frequency j/-mode, with A = —4. For typical ry's of order 1, 
(5i(l) does not vary much from its zero velocity value of 2. The resulting curve A(i!)/Ag starts linearly a,t v = 
from 1 with slope r] and diverges at t; = 1. Thus, in the infinite N limit, the velocity never exceeds the wave-speed. 
At any finite A^, however, the velocity crosses the wave speed at a A of order N'^^^Aq- Since the divergence with A'^ 
is so weak, crossing the wave-speed barrier may not be as difficult as one would naively think. This is especially true 
for small dissipation, where the critical A scales as {rjNY^^. This appears a more likely mechanism for explaining the 
experimental observation of supersonic cracks than the timc-dcpcndcnt forcing hypothesis of Slcpyan. [16] 

The basic structure is unchanged when we go over to the full lattice model. The essential difference is that the lattice 
dispersion relation is nonpolynomial and has an infinite number of positive (real-part) solutions for each eigenmode 
m. The A — v relationship is 

A = VfcFTT TT + ^"Q^-"--; (11) 

Ag (5l,n,m(l + r]vqi,n,m) 

where now the product extends over all positive real-part roots Qi,n,m of the lattice dispersion relation 

= (1 + r?uQ)(4sinh2(Q/2) + A„) - v^Q^ (12) 

for each A„ (Am in the case of gi,n,m)- For a given m, there is one real positive root, Qi,o,m (9i,o,m), and an infinite 
series of complex-conjugate pairs of complex roots, ordered by increasing imaginary part. For large n, the imaginary 
part increases by roughly 27r for each successive root. 

Again, for symmetric cracks we can evaluate analytically the macroscopic (large A'') limit. We obtain 



A = (1 _ ^;2x-l/4 / 2(1 + 7?i;goc,o) 

Ag ^ V ^-^.o 



go,n(l + rivqoo,n) 
9oo,n(l +'r]vqo,n) 



n 



1/2 



(13) 



where qca.n is the root corresponding to the highest frequency A — —4 eigenmode and plays the role of Qi{l) of the 
previous x-continuum result. The go,n are the roots corresponding to the A = eigenmode. These do not have a 
counterpart in the x-continuum calculation as the n = real solution vanishes, and only the lattice-induced n ^ 
modes enter. 

As indicated by the way we expressed this result, we can consider it as essentially the a;-continuum result, Eq. (11), 
with the real lattice g'00,0 replacing Qi(l), modified by a multiplicative correction factor involving the complex lattice 
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FIG. 1. V vs. A/Ag in the x-continuum approximation, Eq. (11), and in the exact lattice model, together with the lattice 
result truncated after the n = term, |n| = 1 and \n\ = 5 terms. 



modes. To understand the usefulness of this way of thinking, as well as its limitations, we present in Fig. 1, for 77 — .5, 
the exact numerically computed relationship Eq. (13), along with the x-continuum result Eq. (11). In addition, we 
plot the lattice result, truncated after its n = 0, |n| = 1 and |n| = 5 terms. We see that at larger velocities all these 
results are close, indicating that the lattice-induced shift in q as well as the additional lattice modes play little role 
at these velocities. At smaller velocities, the various approximations differ significantly from each other and from the 
exact curve. We see, in fact, that as v approaches 0, more and more terms must be included in the product to achieve 
an accurate result. The calculation of the limiting behavior at small velocities requires summing all the terms. The 

result of the calculation is that for all 77, as w — > 0, A approaches A|q+ = \/ 1 + \/2Ag, the maximal A for which an 
arrested crack exists. This generalizes the result of Slepyan for infinitesimal dissipation. As v increases, A decreases 
linearly with the 77-independent slope, — A|q+ /2 so that the bifurcation from the arrested is subcritical and universal. 

More progress can be made in the large 77 limit. Here, at fixed A, the velocity goes to zero as 77 increases, so that 
the ratio (p = rjv is fixed. In this limit, we can calculate the infinite product and find 



A 
A^ 



1 1 

coth( — ) + V2 
Z(p 



(14) 



This infinite-77 result, together with the exact result for various 77's, is presented in Fig. 2. 

We see that this calculation does not reproduce the subcritical bifurcation from the arrested crack at small velocities, 
which is a higher order effect. We can evaluate this I/7; correction near the bifurcation at small 0, and find 
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A.A|.,(l + ^.-'*)(l-|;) (15) 

which reproduces the small v behavior described above and shows that the 77 dependent corrections are in fact 
exponentially small in v. The resulting A — v curve starts at A|o+ at w = 0, heads back linearly for a short distance 
of order 1 777 (In 77)^ and then sharply veers forward. 

A last result worth noting is that whereas the Kelvin viscosity model analyzed herein has a nice macroscopic limit 
when expressed in terms of Aq, the model with Stokes viscosity, where the dissipation in put in the masses and not in 
the bonds, does not have such a limit. There an 0(1) Stokes viscosity at the microscopic level changes the continuum 
elastic fields and requires an ever- increasing A/Aq as the sample is made wider. 



IV. THE x-CONTINUUM MODEL 



We begin our analysis with the solution of the x-continuum model, Eq. (6), introduced in Kessler, et. al. [11]. It 
is important to remember that is this model, the lattice structure in the y-direction is left unchanged. The solution 
of the lattice model is similar in structure to that of the x-continuum model, but the latter is a simpler context in 



4 




FIG. 2. rjv vs. A/ Ac for rj — 2,4:, 8, 16 along with the asymptotic result for large rj, Eq. (14). 

which to develop the necessary techniques. Furthermore, the x-continuum model is an interesting approximation in 
its own right, which captures a significant amount of the structure of the full lattice problem. 

In Kessler, et. al. [11], a Wiener-Hopf analysis of the problem was initiated. In this analysis, the key technique is 
to decompose all the terms in the steady-state equation of motion into terms analytic in the upper- and lower-half 
planes respectively. However, the analysis was not carried to completion, due to the presence of one term whose 
decomposition was not evident. Here we use a trick to accomplish the decomposition of this last remaining term, 
and thereby complete the solution of the problem. We choose not to reproduce the lengthy preliminary stages of this 
calculation, for which the interested reader is referred to [11]. We do however reiterate the definition of the relevant 
notations introduced there, so that the current exposition is minimally self-contained. 

The problematic term, from [11] Eq. (42), is 



ik 



(16) 



where the x, Q, and Q are the roots of a certain family of cubic polynomials. In detail, let £i, {I — 1, . . . , N — 1) he 
the eigenvalues of the {N — 1) x (TV — 1) coupling matrix A^Ar_i(l) defined in Eq. (4) above. Define the polynomial 
specifying the dispersion relation, P{X,Q), by 



P{X, Q) = rivQ"^ + (1 - w2)g2 + (1 + T^yQ)x 



(17) 



Then, P{im,Q) has, for each m, three roots, one positive which we denote Xi,m, and two with negative real parts, 
which we denote by — X2,m, —X3,m, so that all the x's have positive real parts. Similarly, denote the eigenvalues of the 
N X N matrix MN{k) by A™, m = 1, . . . , A^. Then Qi^m, —Q2,m, and — Qs,™ are the roots of P{Am, Q). Likewise, 
denote the eigenvalues of M{0) by A^. Then qi^m, —<l2,m, and —q3,m are the roots of P{Xm, Q)- 

It is apparent from these brief remarks that the troublesome term in its current state involves singularities and poles 
in both the upper- and lower-half planes. To proceed, we rewrite the numerator using the following manipulations: 



'[[{K-ixi,i){K + iX2,i){K + iX3j) 



1 — irjvK 
irjv 

1 / 1 — irjvK 
k 



N-l 



ITjV 

iijv 



1 



fc \ 1 — irjvK 



detN-i{f{K)I + M{l)) 



N-l 



[detN {f{K)I + Mm - detN {f{K)I + M{k))] 
W{K - iqi^,n){K + iq2,m){K + iqz,m) 



5 



1[{K - iQi,m){K + iQ2,m){K + iQs:^ 



(18) 



The first line of this cliain employed an identity from [11], Eq. (40), relating the numerator to the determinant of a 
certain matrix formed from A^jv-i(l) and the identity matrix T together with the function 



f{K) = [ir]vK^ - (1 - v'^)K'^] /(I - ir]vK). 



(19) 



The second line in this chain, claiming this determinant is equivalent, up to a constant factor, to the difference of two 
N X N determinants can be proven by expanding each of the matrices about the first row. The last line reexpresses 
each of these two determinants using more identities from [11], Eqs. (38-39). 
After these manipulations, our term can be written 



UijK - iXi.i){K + iX2,i){K + ixz,i) 

]\^{K - iqi,m)iK + lQ2,rn){K + iQs^m) 



rjv 



1 — ir)vK 



n 



K 



n 



(K + iq2.m){K + ig3,m) 
(i^ + i<92,m)(i^ + i<93,m) 



(20) 



As the outside factor has a pole at —i/rjv in the lower-half plane, the second term only has singularities and poles 
in the lower-half plane and so is in the desired form. The first term is still mixed and requires further massaging. The 
idea is to subtract out the unique lower-half plane pole so that what is left has only upper-half plane poles and zeros. 
Thus, 



K - iQi, 
1 — ir]vK K — iqi^r 



n 



r]v 



1 



imK — + ai 



(21) 



where now g~ has only upper-half plane poles and zeros. We will not need the explicit form of g~ in the calculation. 
What we have is now sufficient to solve for , the Fourier transform of the displacement of the bottom masses in 
the crack region, ui{x)9{x). Using [11], Eq. (42), we find 







{K + ){K + iqs,,n) 



{K + lQ2,rn)iK + lQ-S, 



K 



n 



+«i(0) 



1 — ir}vK 



1 -I- -qvQi 



*0+ Q2,mQ3,m 

{K + iq2,^){K + 

^q3,m) 



n 



{K + iQ2,m){K + iQ3,m) 



(22) 



Solving for u~^, we find 



u+ = 



K + iO+ Q 



n 



{K + iQ2.,m){K + iQ^^rn) 



{K + ){K + 

(1 -I- 'qvQijn){K + lQ2.,m){K + iQ^Au) 



TT — 

(1 



{l+rivqi^m)iK + )iK + iq3^rn) 
Fourier transforming and evaluating at x = 0+ yields 

Wl(0) = A[[ 

SO that 



(23) 



Q2,mQ3,- 



-ui(0) 



"L {^ + Vvqi,m) 



^.i(o)=An 



q2.n,M:',j„('i- + >lt''ll.m) 
(32,mQ3,m(l + VvQl,m) 



(24) 



(25) 



Using ui(0) = e, Ag = e\/kN + 1 and the relations (see [11] Eq. (43)) 
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m 



we obtain our desired result 



(26a) 

(26b) 

(27) 



The primary benefit of this method of solution over the direct approach employed in [11] is that for the symmetric 
crack [k = 2) we can take the large- limit. To do this, we break up the iV— fold product into two terms. The first is 



(28) 



We transform the product into the exponential of a sum over logarithms, a sum which for N large we can approximate 
by an integral, via the Euler-MacLauren Summation Formula (EMSF) [14]. Thus 

r-N 

lnni«/ dm[\u{l + rivQi^m)-M'^ + nvqi,m)] (29) 
Jo 

Now, for k = 2, = —Asiia^^ ^N+i ) ^^nd = — 4sin^( '^^^^^^'' ). If we define a = m/N, then we see that 



gi(a) = Qi(a-l/(2iV))«Qi(a) 



1 dQi 
2N da 



The integral is now a total derivative, and so 



InHi « i / da^ [ln(l + vvQi,m)] = J [ln(l + VvQi{l)) - ln(l + VvQi 

2 Jq da 2 



(0))] 



(30) 



(31) 



When a = 1, m = N and so w —4, and so Qi(l) satisfies 

= P(-4, Oi(l)) = vvQiilf + (1 - v')Oi(l)' - 4(1 + VvQiil)) (32) 

Similarly, when a approaches 0, so does m and so also A„. This in turn implies that Qi{0) = 0. So, finally, 

ni«[l + r?t;Qi(l)]V2 (33) 

The second factor is slightly more difficult to treat, since the numerator and denominator both vanish as m — > 0. 
To handle this, we regularize the product by multiplying and dividing by Ylm v^Vra/A^, which we can perform 
analytically. Then, the regularized product, n|^, can be transformed to the exponential of an integral of a total 
derivative, which can be calculated explicitly. In detail. 



Ql,m,'\/~^v 
9l,m\/ — Arr 



(34) 



so that 



In 



In- 



2 da ^-k{a) 
In (Qi(l)/2)-ln (l/^/l^^ 



where we have used the fact that for a small, Qi{a) k, ^ — A(a)/(1 — v^). Thus 

1/2 



Qi(i)vT^ 



(35) 



(36) 
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Also, 



U2 = U^l[^/Am/Xm 

m 

= U^^/detM{2)/detM{0) 
= UfV2N + 1 



so putting all the pieces together yields the simple result 



.ni 



/ 2(l + r?,;Qi(l)) 



(37) 



(38) 



which expresses A in terms of Qi(l), the wave-vector at the end of the Brillioun zone. 

The most striking lesson of this formula is that A diverges at u = 1, the wave speed. Thus, while at any finite A'', 
there is no upper limit to the velocity, at infinite N the wave speed is an absolute upper bound to the crack velocity. 
At large A, v approaches unity from below as 1/A^. A second lesson is that at small velocity, A Ag{1 + r^v), so 
that A approaches Ag linearly, as is generally true for this x-continuum model. A third implication is the behavior 
at large r}. For fixed A, v decreases as r] gets large, so that Qi{l) satisfies 



« nvQi{if + Qi{if - 4(1 + vvQi{i)) = (1 + vvQiimQiiif - 4) 



so that (3i(l) ~ 2. Substituting this in Eq. (38) gives 



(39) 



(40) 



or 



(41) 



In this large 77 limit, of course, A is a fimction of the scaling variable rjv, which was first introduced in [15]. 

We have seen how at infinite iV, the crack speed v never crosses unity, the wave speed. However, at any finite 
N, there is a A for which the crack speed crosses unity, which must diverge with N. We now calculate how this 
threshold scales with A''. The key to the calculation is 112, since it is the vanishing of 112 which leads to the divergence 
of A at w = 1 for AT infinite. To compute the value of 112 at = 1 for finite large N, we need to choose a difii'erent 
regularization. We now define 



n 



9l,m(-A,„)V3 



(42) 



so that 



Ho = m 



n 



A. 



1/3 



= {2N + iy/^ii. 



(43) 



Now, since (5i(q;)/(— A(a))^/'^ approaches the finite limit l/rj^^^ as a goes to 0, 11^ has a finite limit at i; = 1 as A'' 
goes to infinity, namely II2 ~ ^/Ql{l){r]/Ay/^. Using the infinite N limit of Hi from Eq. (33), we get 



A = V2NTT^ ~ (AT) 
Ag II2 



1/6 



2(l + r,0i(l)) 



Qi(l) 



1/2 



(44) 



Thus, the threshold A scales as N^^^Aq, in accord with the numerical evidence discussed in [11]. The coefficient goes 
to 2 for T] large, and vanishes as 77^^^ for small rj. 

Another manifestation of this same phenomenon, the disappearance of the v = 1 crossing in the infinite- AT limit, is 
the nonuniformity of the the large- A'' limit as v approaches 1. Working out the corrections to the EMSF, we find that 



ni(Af) «ni(Ar = oo)* 1- 



■nrjv 
8\/l - v'^N 



(45) 



8 



and 



n«(iv 



oo) * I 1 + 



rjv 



16(1 - u2)3/2A/- 



(46) 



Thus, the relative error of the infinite-A^ approximation is 0{1/N), and diverges as v approaches 1 as (1 — u)^'^/^. 
It is also interesting to note that the relative error vanishes as v goes to zero, son that the infinite-iV approximation 
becomes better at small velocities. 

One last interesting piece of information we can derive from our solution is the size of the "process zone" , the region 
where the solution from continuum elastic theory breaks down. The leading-order macroscopic solution was derived 
in [11], and exhibited the classic square- root singularity at the crack tip, x = 0. This singularity in really present only 
at infinite A'', and is cut off by the upper limit on the Q's, (relative to the smallest Q ~ 1/-^) at finite A^. We can 
determine the structure of the process zone which replaces the singularity by studying our exact solution for w+, Eq. 
(23), for i^'s of order 1. Using the EMSF to evalute the infinite product, similar to the derivations above, we find 



iA 



Qiil)VT^{K + iQ2{mK + tQ3{l)) 



1/2 



y/2N+l{K + iO+) 



+ 



le 



K + i/riv 



2{K + iQ2{0)){K + iQsiO)) 

(1 + rivQi{l)){K + iQ2{l)){K + iQsjl)) 
(1 + r}vQi{0)){K + iQ2{0)){K + iQaiO)) 



1/2 



Using Qi{0) = 02(0) = 0, Qs{0) = (1 - v^)/r]v, and the result for A, Eq. (38), we get 



(1 + r)vQiil))iK + iQ2imK + iQ^jl)) 
{K + iO+){K + i{l-v'^)/rjv) 



1/2 



1 



1 



iO+ K + i/r]vJ K + i/ijv 



(47) 



(48) 



Examining this expression for small K, we find the expected K~^^^ singularity, which gives rise to the square-root 
singularity of the outer solution. The coefficient of the K~^/'^ singularity to leading order in N is Ai'^/'^N~^I'^{\ — 
^2^-1/4^ which reproduces the same 77 independent coefficient of the square-root singularity, or equivalently stress- 
intensity factor, found in [11]. The structure of the process zone is governed by the other singularities in ■u"*' that lie 
off the origin. In particular, the size of the process zone is determined by the singularity nearest the real line. For 
small ?7, this is at K = — iQ2(l) ~ — 2?7\/l — v'^ : so the process zone is truly microscopic, unless the velocity is very 
close to 1. For large 77, the dominant singularity is at iiT = —iQai^) ~ —i/rjv, so the process zone grows linearly with 
T] in size. 



V. THE LATTICE MODEL 



In this section, we generalize our solution of the continuum model to the lattice model. For ease of presentation, 
we will present the derivation only in the iV = 1 case. The case of general follows in a straightforward manner 
from this derivation and that of the continuum finite A^ model presented in the previous section. 

Our derivation follows directly along the lines of our WH treatment of the continuum N = 1 problem in [11]. The 
equation of motion of the steady-state crack is 

uit) = (1 + 7?^) [uit + 1/v) - Suit) + u{t - l/v)] - ke{-t){l + Tl^)u{t) (49) 
Upon Fourier transforming, we find 

= (1 - i-nvK) ^4sinh^(^) - 1^ w + v'^K'^u - k{l - ir]vK)u- + A6{K) - kr]vu{0) (50) 

where u is the Fourier transform of u and are the transforms of 6{±t)u{t). We define the function 

R{X; Q) = (1 + r]vQ){'ismh^{Q/2) + A) - v^Q^ (51) 

in terms of which 

= R{-{1 + k); -iK)u- + R{-1; -iK)u+ + AS{K) - kr]vu{0) (52) 
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This function, R{\: Q), which is the lattice equivalent of the polynomial P employed in the previous section, has 
not 3 roots, but in fact an infinite set of zeros in the complex plane. We shall label these zeros according to their 
real parts, (5i,n (9i,n) are the zeros of — (1 + k);Q) 1;Q)) with positive real parts, and Q2,n' {<l2,n') are their 
counterparts with negative real parts. The indices n, n' run over the entire infinite set of zeros but are otherwise left 
unspecified for now. We can decompose R in terms of its zeros 



R{-{l + k);-iK) 
R{-l;-iK) 



K 

Ql,n 



)(1 



I- — ) 



Q2,n 



qi,n Q2,n' 



(53) 



Using this, we rewrite the equation of motion: 

= -(l + A:)n 



qi,n{K - i<3l,n) ~- 

u 



Qi,n{K-iqi,n 



TT Q2,n'(.K + iq2,n') ,_)_ 
q2,n'{K + iQ2,n')^ 



(54) 



As in the last section, the hard part is to decompose the last term. The trick is the same, rewriting the numerator 
as the difference of i?'s. 



'Ha 



1 



1 R{-l;-iK)~ R{-{l + k): -iK) 



K 
' 51,, 



_ 1 + fc -J-j- gi,n(i^ - iQl.n) 1 



Q2,n'{K + Zg2,n') 

1 - i'qvK J-J- Qi.n(,K - iqi,^) 1 - i'qvK J-f q2.n'{K + iQ2,n') 



n 



(55) 



As before, the second term is now fine, but the first term is still mixed. Again we subtract out the unique pole in the 
lower-half plane which is what we need to find w+. 



ii ' 



1 



1 - irjvK i-i- Qi,n{K - iqi,n) 1 - ivvK if Qi,„(;i; + %,„) 



n 



qi,n{:^ + Ql,n) 



+ 9~ 



(56) 



where g only has poles and zeros in the upper-half plane. Separating out the pieces analytic in the upper-half plane 
yields 



: 



n 



Q2,n'{K + iq2,n') - + 

u 

{K + iQ2,n') 



+ 



K + iO+ 



— rivu{0) 



1 + fc qi,n{^ + Ql,n) 1 -TT (52,n'(-fl' + ^92,71' 

1 - irjvK J-J- Qi,„(-^ + iqi,„) 1 - ivvK -If q2,n'{K + iQ2,n' 



(57) 



Solving for u+ yields 



u+ 



TT g2,n'(-ft' + ^Q2,n') 
K + iQ+^Q2,n'{K + iq2,n') 

rjv 



-wi(0) 



1 — irjvK 



{l + k)Yl + ■nvQi.n){K + iQ2,n') _ ^ 



, Ql,nQ2,n{^ + ){K + iq2,n') 



Fourier transforming and evaluating at a; = 0+, we find 



(58) 



ni(0) = An 



92, ra' 
f Q2,n' 



Ul{0) 



(1 + ^)11 



(1 + rjvQi.n) 

-, Ql,nQ2,n'{^ + VVqi,m) 



(59) 



so that 
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Vvqi,n) 



(60) 



As A(3 = ui(0)\/l + k, we obtain our desired result 



jl+rjvQi^n) 

Ql,n{^ + VVqi,n) 



(61) 



As there is exactly one real positive root of i?(A; Q), it is convenient to assign this the index and to label the complex 
roots in order of imaginary part, so that for example and Qi^-n are complex conjugates. It is clear the basic 
structure of the lattice result is similar to the continuum result Eq. (27) above, with the continuum Qi^i, replaced 
by their lattice counterparts Qi^o, <Zi,Oi £^nd multiplied by a correction factor due to the additional infinite hierarchy 
of complex Q, g's which solve the lattice dispersion relation. 

The generalization to finite A'' is straightforward and is left as an exercise to the reader. The result is the direct 
generalization of the = 1 result. At finite A^, there is a set of zeros with positive real part of R{Am', Q), {R{^m.', Q))i 
for each m = 1, . . . , A^, now labeled Qi,n,m {Qi,n,m)- Then 



— = VkN + Tl[ 



.n.rn ) 



(62) 



is the solution to the lattice problem at finite N. It of course reduces in the limit — > 0+ to the result of Marder and 
Gross [9]. 

As in the continuum, this rather unwieldy formula simplifies tremendously in the symmetric crack case fc = 2 as 

A^ goes to infinity. The procedure for evaluating the limit is similar to the continuum calculation and so wc do not 
present the details. What enters again are Qi,„(a), at the two extremes of the Brillouin zone a = 0, 1. If we label 
Qi,n(l) = <loo,n, <3i,n(0) = qo,n then they satisfy the dispersion relations 



= i?(-4;goo,«) = (1 + r7W9oo,«)(4sinh^((7oo,«/2) -4) -w^gr, 
O = i?(0;go,«) = (l+r7^'9o,n)(4sinh2(go,„/2))- 1.2^2^ 



2.2 
CJO.n 



In terms of these g's, the infinite A^ limit solution is 



A 



(l_^2)-i/4 /2(T+,ay.^,o) 

y 9oo,o 



n 



(1 + vmo,n) 



1/2 



(63a) 
(63b) 



(64) 



Again, this is very essentially similar to its continuum counterpart, with the real lattice wave-vector goo,o playing the 
role of the continuum wave vector Qi{l), and with a multiplicative correction due to the presence of complex lattice 
wave- vectors. It should also be noted that this result reduces to that of Slepyan [6] in the ^ 0+ limit. 



VI. THE SMALL VELOCITY LIMIT 



Wc begin our explorations of the content of our key result, Eq. (64) by examining the t] fixed, u ^ 0"*" limit. It is 
not sufficient to simply set v = 0, since as v gets smaller, more and more terms contribute significantly to the infinite 
product, as seen in Fig. 1. The proper treatment is to replace the infinite product by an infinite sum of logarithms 
and then approximate the infinite sum by an integral via the EMSF. Note that for v = 0, qoo,n satisfies sinh^ _ 
with the solution 

q^^°=27Tin + u, (65a) 
where w is the unique real root of the equation, namely co = 21n(l -|- \/2). Similarly, 

q^oTn = 27rm. (65b) 

As we discussed above, we need to consider w — > 0, n — > oo, a = 2'jTr]vn fixed. Then, writing „ = 2nin + ujao, ijJoo 
satisfies 
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. 1 2 '^oo , 

smn — — = 1 



(27rm + bj^oYv'^ 
4(1 + r]v{2nin + lOoc)) 



477^(1 + ia) ' 



Similarly, 



a 



2 4r?2(i + ia)' 
We can now easily approximate the first infinite product, 

TT 77^ + Wqoo,n 

^1= 11 Y 



yielding 



In 



da 
2'Kr}v 



[ln(l +ia + r]vw^) — ln(l + + rjvujQ)] 



(66) 
(67) 

(68) 

(69) 



7-00 277 1 



■ Wo 



(70) 



The second product is somewhat trickier, because a naive expansion diverges at small a. We define a regularized 
product 



(7oc^„(27rw) 
go,n(27rm + w)' 



(71) 



so that (using the product formula for sink, and the fact that sinh(a;/2) = 1) 

2nin + u) 



n2 = n«n 



n#0 



27rm 



^ sinh(a;/2) 

w/2 2 

= -nf 



(72) 



Our regularized product is now easily approximated, 



In 



f da r / 
nfw / In 



ia + WocV^ 
ia + ojriv 



In 



ia + LUorjv 



la 



f 

J — c 



da Woo — w — Wo 
27r ia 



Using the identity 



we obtain our desired result 



da- — = TTW, 

, 1 +ta 



(73) 



(74) 



A 



0+ 



2ni 



-,1/2 



9oo,0n2 



e*^/* exp [ - 



da Woo — w — Wo 
27r a(l+ia) 



(75) 
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This result is, as desired, explicitly independent of but would appear to depend on rj through the very nontrivial 
7] dependence of Wqo and ojq under the integral. It is possible to explicitly evaluate the integral for small and large rj. 
For large r], Wqo — w 0{l/rf) and wq ^ 0{l/ri), so the integral vanishes and so 



A|o+/AG = e'"/4 = VI + 1-554. (76) 

For small r] [6], cOoo — is concentrated at small a ^ 0{r]), so it is appropriate to convert the integral into a principal 
value integral and do the u) integral immediately. In the remaining integral, we change variables to /3 = a/2ry. Then, 
the denominator in the integrand reduces to 1//3, so only the odd (i.e. imaginary) part of Woo — <^o contributes. For 
/? > we find 

f3 < 1 

Imuoc = 2Im sinh"i(v'l - P"^) = { 2sm-^ y/ - 1 1 < /? < 2 (77) 



and 



The integral thus becomes 



/mco = 2/m sinh-i(^) = | ^ f { 



(78) 



r°° da oj^ - w - ujQ _ ioj f^dp sin-^p /"^ d/? sin" V/j^ - 1 - 7r/2 
J_^2n a{l + ia) ~ 2 Vo ^ ^ Vi /3 

= i ln(l + x/2) - i In 2 + I In(^^t^) 

= (79) 

So, again the integral vanishes, and the result for large and small r] is the same. One is lead to guess that in fact 
the integral vanishes for all rj, as indeed a numerical computation confirms. This is physically reasonable, since A[o+ 
should be nothing other than the maximal A for an arrested crack, which was previously found numerically [11] to be 
approximately 1.55Ag. This maximal arrested crack A is the result of a static calculation, and is of course completely 
independent of r]. The vanishing of the integral can be demonstrated analytically and is the result of the fact that the 
integral has no singularities in the lower-half-plane. One can then close the contour there and the result is identically 
zero. 

To see this, one has to study the analtyic structure of the functions (jjoo(a), uJo{a). Consider first uJo{a) = 
2sinh~^(2/(a)), where y^{a) = — a^/(477^(l + ia)). Since sinh~^{y) = 2ln{y/{y'^) + yTTy^, loq has a branch cut 
singularity along the line y'^ = —r where 1 > r > 0. Working out the algebra, in the complex a plane this works out to 

be, for ?7 > 1, two separate curves. The first is a segment along the upper imaginary axis from a ~ 2ii](r] — \/rp~-i) 

up to a = 2i'q{r]+ sjrf — 1). The second is the circle of radius 1 centered at the point r] = i. Similarly, Woo has branch 
cuts for 2 > r > 1, which are two finite segments along the positive imaginary axis extending above and below the wq 
branch cuts. For y/2/2rj < 1, the branch cut for cjq is a sector of the circle, while the branch cut for Woo is the rest of 
the circle and a finite piece of the the entire imaginary axis extending centered about 2i. For 77 < V2/2, the branch 
cuts are confined entirely to a part of the circle. Thus, the singularities for all 77 lie entirely in the upper-half plane 
and, as advertised, the integrand is analytic in the lower-half-plane and so the integral vanishes. 

The next step is to extend this calculation to next order in v. There are two sources for this first correction in 
V. One comes from the higher-order velocity dependence of qoo,n, qo,n- The other comes from the EMSF correction 
to the replacement of the infinite sum by an integral. The calculation of the first piece is similar in structure to the 
leading order calculation, just more involved. We expand qs,n ~ 27rin + lOs + vag, where the subscript s = 00, 0, and 
find 



auJs 2i-a 

= ^ — — 1 , ■ 2 - ^0 
2r]smhu)s 1 + »Q: 



It is straightforward, though tedious, to substitute this in Hi , 11^ and expand, giving a multiplicative correction factor 
of 



iv da 
^ 2" / ^ 2^ 



(81) 
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FIG. 3. V vs. A/Ag for 77 = 0.5, 1, 2 along with the asymptotic resuh for smaU v, Eq. (83). 

The integrand is again a very nontrivial function of 77 and a, but again a miracle occurs and the integral vanishes 
identically (as seen by numerical computation) for all 77! The same analytic argument as above can be used to prove 
this point. 

This leaves us with only the second source for a 0{v) correction, namely the first endpoint EMSF correction to the 
integral. There are only exponentially small corrections to the integral representation of IIi , but since II2 does not 
include an n = term, we need to subtract the n = limit of the summand, namely 

27rin , , 1 , , , 

lim In ) = In « -v, 82) 

from Innf^. This gives a multiplicative correction factor of (1 + v) to Xlf , so we find that for small velocity 

A|o+(l-t;/2) (83) 

Thus the leading small v behavior of A is completely independent of rj. However, it has A as a strictly decreasing 
function of v. As we shall see in the next section, the turnaround for larger v is a nonperturbative effect. For now, we 
will conclude this section by showing in Fig. 3 a plot of the small velocity region of the graph for various 77's, together 
with our analytic approximation. We see that the analytic result is confirmed. 



VII. LARGE rj LIMIT 

We now turn to a study of the large limit. In this limit, as first pointed out by Pla, et. al. [15], the relevant 
variable is rjv. Thus, we study the limit rj — s- 00, v ^ 0, (j) = rjv fixed. As we shall see, this calculation will shed much 
light on the small v results we obtained in the previous section. 

To begin the calculation, we need the goo.n's and go.n's at t; = that we obtained in the previous section, Eq. (65). 
Then 

^ -^i l + r]vqo,n \ + (t>{2mn) sinh(^) ^ ' 



n— — 00 



Similarly, 

900 n TT 27rm + UJ 2 . UJ 2 
112 = - " " 

Thus, 



^ T-r goo.n 77 imn + uj I . UJ I 

Ha = I I W — — = - smh(-) = - (85) 
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A 


2ni 


1/2 


rsinh(l±|^)l 


1/2 




.9cx>,0n2. 




. sinh(^) 





coth(-i-) + V2 



1/2 



(86) 



We can invert this relation, solving for (f) in terms of A, which yields 



in 

, (A) 



(87) 



For large A, this approaches ^ [(A/Ag)^ — \/2]. This asymptotic result, which is also presented in Fig. 2, is 
to be contrasted with the result of our continuum calculation, where wc found (/) = i [(A/Ac)^ — l] . Thns, the 
continuum infinite-?? calculation for all A essentially reproduces the large-A limit of the lattice calculation, with the 
correct functional dependence, but with the graph just shifted down slightly. It is also worth noting that including 
just the n = term, instead of the whole infinite product, also gives the same result, with an intercept of 2/uj which 
is intermediate between the continuum calculation and the exact asymptotic result. As A decreases, the true r] = oo 
curve falls below the asymptotic result, so as to intercept the A-axis at A|o+. The approach is singular, as can be 
seen by looking at Eq. (86) for small a. We find 



A~ A|o+(l + 



A2 

A|2 



(88) 



with an essential singularity at small (j). 

Examining Fig. 2 more carefully, we see that our infinite rj result has failed to capture one of the most salient 
features of the finite ry data, namely the subcritical nature of the bifurcation from the arrested state. Instead, it 
possesses a (very-) marginally supercritical onset of the moving crack. To reproduce the subcritical bifurcation from 
our analytics, we need to generate the next order correction in I/tj. 

We begin by generating the next order correction to the g's. We find that goo.n does not change to this order, but 
now go n ~ 2mn + where 



Wo 



2TTin(j) 



- 1 tan ^ 27rn<^ 



77(l + 47r2n2^2)i/4 
This induces a multiplicative correction to A of 



1-E 



n#0 



47rin^(l + 27rm(/)) 



(89) 



(90) 



We are interested in the effect of this correction at small cj), in which case we are again free to replace the sum by an 
integral. If we add in the n = term to the sum, the error will be exponentially small in l/cp. So, up to exponentially 
small terms, the correction for small <p is (defining a = 27rn^) 



da 



1 



27rr?(l + a2)V4(i+iQ,) 



(91) 



The integral vanishes, as can be seen by a substitution of variables x = (1 + a^)-!/^^ f^ct, the integral is nothing 
more than the first-order expansion in l/ry of the integral in Eq. (75) which we found vanishes identically in rj. We are 
thus left with a correction factor of simply (1 — (l)/2r]) = (1 — v/2) up to exponentially small terms. This is precisely 
the small v correction we found in the previous section. The full behavior to this order for small is thus 



A|o+(l- 



A^' 
^lo+ 



,-1/0 



)(1--) 



(92) 



This has the subcritical bifurcation we are seeking. As increases from 0, A decreases from A|o+ due to the influence 
of the second factor, until the exponential kicks in and causes A to turn around and start increasing. The 4> at which 
the turn-around occurs is, for large rj, of order 1/lnr? (translating to a velocity of order l/rjlnry) which goes to as 
r] goes to oo, but very slowly. Thus at infinite 77 there is no turnaround and A strictly increases with (f> as we found 
in the zeroth-order calculation at the beginning of this section. The minimum A lies, for large rj, an amount of order 
l/?7(lnj7)2 below A|o+. 
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Thus wc sec that it is the subdominant pieces that are responsible for the increase of A with v, while the perturbative 
pieces give rise to the subcritical bifurcation. Analyzing the subdominant pieces in a little more depth, it is easy to 
see that for r] > \/2/2 the leading subdominant piece goes as exp(— 1/771,'). For smaller -q, the subdominant piece falls 
less rapidly, and has an oscillating component, due to the off-axis branch cut assuming dominance. This picture is 
consistent with the numerical evidence. 



VIII. STOKES VISCOSITY 



It is worthwhile to contrast the behavior we have seen for Kelvin viscosity with that which obtains for Stokes 
viscosity, where the dissipation is associated with the mass points and not the bonds. The calculation in this case 

is much simpler, since the troublesome r] term is not present. For our purposes, it is sufficient to consider our 
a;-continuum theory, as the conclusions we obtain carry over to the full lattice model. The result for m+ is 



iA -r-r (72,m(^ + «(32,m) 



where now the Q's satisfy the dispersion relation 

(l-u2)Q^-6wQ„ + A„ = (94) 

(and the q's the parallel form with A^) and b is the Stokes viscosity. This can be seen by a simple limiting procedure 

applied to Eq. (23), or by replaying the derivation leading up to Eq. (41) of [11] with b instead of 77. This form of the 
solution can be shown to be equivalent to that obtained by Marder and Gross [9]. This result leads to the solution 
for A: 

A = eTT^ (95) 



We are interested in the large- A/' limit, which we obtain be defining the renormalized product 

n« = TT Q^''"/-^) (96) 

\^ ?2,m(-A„) 

since the Q2's are linear in A for small A. Applying the EMSF, we find that for large A'', 



so that 

A^{2N+ l)e{b\^ + 16(1 - i;2))V4yC^ (98) 

The key difference between this formula and the parallel one for 77 is that A/e is proportional to N, and not N^^"^ 
as before. The reason for this is that the Stokes viscosity is most effective at damping small wavelengths, and so 
affects the macroscopic stress fields. The Kelvin viscosity does not damp out small wavelengths and only acts on 
short wavelengths. Another way to see this is to compute the stress intensity factor, which in the Stokes case is 
inversely proportional to Vb. The driving force required to propagate the crack is thus much larger in the Stokes 
case. In particular, in the Stokes case there is no macroscopic scaling limit, where things just scale with the Griffith 
driving, Aq. For these reasons, we feel that the Stokes viscosity is not a good model of dissipation for studying crack 
propagation. 

The only way to obtain a nice macroscopic limit where A scales like Aq is to artificially scale b with N so that 
h = bf)/N. However, this procedure has no physically satisifying motivation, especially when the Kelvin viscosity 
model suffers none of these defects. 
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IX. CONCLUDING REMARKS 



We close by making a few comments about this work and prospects for future extensions. First it is important to 

note that the present work is Umitcd to a consideration of the steady-state crack. Thus, aside from general issues of 
the size of the process zone, the major output of this problem is the velocity-driving relation. Here the most striking 
qualitative effect of Kelvin viscosity is near threshold, reducing the extent of the backward bifurcation. Significantly 
above threshold, the major role of viscosity is to provide a velocity scale, so that the crack velocity becomes inversely 
proportional to the viscosity. It is important to understand how viscosity impacts on the stability of the crack. It is 
clear, as Marder and Gross have pointed out [9], that the steady-state crack is unstable in the regime of the backward 
bifurcation. The more interesting question is in the higher-velocity regime. Here, no systematic studies have been 
done to examinne the role of viscosity. It is not clear that the piecewise-linear model considered here is altogether 
appropriate for studies of stability, as instabilities can be masked by inconsistencies of the steady-state solution. We 
look forward to reporting on work in this direction soon, along with generalization to the problem of Mode I cracking. 
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